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CT  (Maximum  xuu  wo 


We  develop  an  integrated  approach  to  estimating  velocities  and  azimuths  from  a  collection  of  local  arrays 
and  then  fusing  the  data  into  Bayesian  locations  and  their  associated  uncertainty  ellipses.  A  small-array  theory 
is  developed  that  characterizes  the  performance  of  local  optimal  detectors  under  signal  correlation  and  decor¬ 
relation  scenarios.  We  compare  the  performance  of  maximum  likelihood  estimators  such  as  the  beam  power 
and  the  generalized  beam  power  as  a  function  of  array  geometry  and  Bignal  to  noise  ratio.  Optimal  local-array 
geometries  are  suggested  that  are  relevant  to  the  problem  of  designing  an  optimal  infrasound  array. 

Wave-number  estimators  along  with  estimated  variance  covariance  matrices  are  used  as  input  to  study  the 
size  and  orientation  of  90%  posterior  probability  ellipses  for  various  likely  subsets  of  detecting  stations  within 
the  global  infrasound  array  proposed  for  the  Prototype  International  Data  Center  (PIDC).  Adding  detecting 
stations  decreased  the  size  of  the  90%  ellipse  by  about  10-20%  per  added  station,  whereas  increasing  the  signal  to 
noise  ratio  from  2  to  3  decreased  the  size  of  the  ellipse  by  30-40%.  Adding  an  inner  triangle  to  the  conventional 
1  km  triangular  array  gave  more  modest  reductions  of  7-10%. 
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3  Estimated  large  sample  standard  deviation  of  azimuth  estimators  as  a  func-  9 
tion  of  baseline  distance  using  the  Mack-Flinn  theoretical  coherence  model  for 
decorrelation  with  distance  with  the  solid  lines  for  the  weighted  beamforming 
estimator  and  the  dotted  lines  denoting  simple  beamforming. 

4  Relief  (left  column)  and  contour  (right  column)  plots  of  F-statistics  for  detect-  14 
ing  infrasonic  signal  at  two  bandwidths;  narrow  band  F  shows  peak  at  227(10) 
degrees  and  velocity  c  =  .24(.04)  km/sec. 

5  Relief  (left)  and  contour  (right)  plots  of  beam  power  and  maximum  likelihood  15 
correlated  beam  power  detectors.  Beam  power  shows  peak  at  225(7)  degrees 

and  velocity  V  =  .26(.03)  km/sec.,  correlated  detector  shows  peaks  at  225(6) 
degrees  and  .26(.03)  km/sec. 

6  Possible  prior  distributions  for  standard  deviations  of  measured  wave-number  18 
estimates. 

7  Posterior  probability  (90  and  95%)  ellipses  for  various  array  detection  config-  21 
urations  assuming  7-element  arrays,  perfect  signal  correlation  and  S/N  =  2. 


•  M 
111 


Section  1 
Introduction 


The  monitoring  of  a  comprehensive  test  ban  treaty  (CTBT)  depends  rather  critically  on  data 
collected  from  various  networks  of  recording  devices.  Generally,  there  will  be  collections  of 
small  local  seismic  or  acoustic  arrays,  providing  as  input,  signals  from  a  number  of  sensors 
recording  a  common  event.  For  example,  Figure  1  shows  time  series  from  a  small  array 
of  microbarographic  instruments  that  recorded  a  nuclear  explosion,  detonated  about  25  km 
south  of  Christmas  Island.  The  dimensions  and  configuration  of  the  array  are  unknown  so 
the  pseudo  data  are  constructed  by  aligning  to  a  given  velocity  and  azimuth  assuming  1  km 
sides  on  the  array.  The  problems  of  interest  for  this  local  array  involve  first  detecting  the 
signal  and  then  identifying  parameters  or  features  that  can  be  used  to  locate  the  source  of 
the  event.  It  would  also  be  critical  to  determine  whether  the  detected  signal  belongs  to  a 
benign  class  of  events  such  as  might  be  produced  by  earthquakes  or  mining  explosions  or 
whether  it  might  possibly  be  an  explosion  of  interest  to  the  nuclear  monitoring  community. 
While  the  second  discrimination  aspect  of  the  above  problem  is  of  great  interest,  it  is  the 
first  question,  namely  the  local  array  characteristics  that  influence  location  capability,  that 
has  been  the  focus  of  this  present  contract. 

x  t0*  S«nsor  A 


Figure  1.  Small  Fry  series  after  aligning  and  delaying  to  correspond  with  an  approximate 
velocity  of  .3  km/sec  and  an  azimuth  of  225  degrees  sampled  at  .1  sec. 


If  it  is  supposed  that  a  detection  is  available  by  one  of  the  optimum  processors  to  be  consid- 
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ered  later  in  this  report,  it  will  be  essential  to  have  estimators  for  parameters  related  to  the 
location  of  the  event.  It  is  clear  that  the  angle  of  approach  or  azimuth  and  the  velocity  of  the 
signal  will  be  related  to  its  origin.  In  general,  the  local  array  provides  estimators  for  velocity 
and  azimuth  or  for  two  wave-numbers  0  =  (01,0a)'  that  are  functionally  related  to  both  the 
velocity  and  azimuth,  as  well  as  the  location  x  =  (xi,xa)'  of  the  event.  We  indicate  this 
dependence  by  writing  0\(x),6a(x), . . .  ,0n(x)  for  the  wave-number  vectors  as  they  might  be 
observed  on  n  local  arrays.  The  global  array  consists  of  n  local  arrays  and  the  combining  or 
fusing  the  local  array  estimators  into  an  overall  location  estimator  and  its  uncertainty  is  the 
objective  of  the  analysis. 

As  an  example  of  a  collection  of  small  arrays  monitoring  events,  consider  Figure  2  which 
shows  14  local  arrays  proposed  for  the  southwest  quadrant  of  the  world-wide  network  of 
infrasound  detectors  to  be  incorporated  into  the  International  Monitoring  System,  denoted 
by  IMS  in  the  sequel.  For  orientation  purposes,  the  local  arrays  labeled  1,2, 3, 4, 6  and  11  are 
all  located  on  the  South  American  continent.  Other  stations  are  in  the  oceans  to  the  west 
and  east,  but  still  south  of  the  equator.  The  two  hypothetical  events,  marked  by  +,  are  in 
the  oceans.  Location  ellipses  for  these  events  developed  in  this  report  are  functions  of  the 
geometry  of  the  local  array,  local  signal  to  noise  ratios  for  the  event,  signal  decorrelation 
and  the  time-bandwidth  product.  An  additional  important  factor  influencing  the  size  of  the 
confidence  ellipses  will  be  the  configuration  of  detecting  stations. 
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Figure  2.  Proposed  IMS  infrasound  arrays  in  the  southwest  quadrant  and  two  hypothet¬ 
ical  event  locations  marked  by  +. 


As  indicated  above,  the  objective  of  this  project  is  to  develop  a  theory  that  relates  the 
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local  array  performance  to  the  location  capabilities  of  a  global  network  consisting  of  some 
configuration  of  detecting  local  arrays.  This  requires  a  complete  investigation  of  local  array 
detection  and  estimation  using  maximum  likelihood;  the  details  leading  to  the  equations  in 
Section  2  can  be  found  in  Shumway  et  al  (1999). 

We  begin  in  Section  2  by  estimating  the  wave-number  parameters  under  the  assumption  that 
the  signal  and  noise  spectral  matrices  are  fixed.  As  an  adjunct  to  the  estimation  results, 
we  obtain  maximum  likelihood  estimators  for  the  signal  and  noise  spectra  in  the  perfectly 
corrrelated  case  using  a  frequency  domain  version  of  analysis  of  variance.  In  the  decorrelated 
signal  case,  we  develop  a  version  of  the  EM  algorithm  that  produces  separate  estimators  of 
the  signal  spectral  matrix  and  the  noise  spectrum.  This  will  provide  estimated  values  for 
the  signal  coherence  in  the  decorrelated  case.  As  inputs  for  the  location  procedure,  we  also 
require  estimated  values  for  the  wave-number  parameters,  say  0*  for  the  k  =  1, . . . ,  n  local 
arrays,  along  with  their  estimated  variances  and  covariances.  These  estimated  parameters 
lead  to  maximum  likelihood  estimators  for  velocity  and  azimuth.  The  maximum  likelihood 
estimators  for  the  wave-number  parameters  can  be  obtained  as  the  maximizers  of  either  the 
beam  power  or  the  F-statistic  and  we  give  the  variances  and  covariances  of  these  estima¬ 
tors  under  both  the  perfect  correlation  signal  model  and  for  the  decorrelated  signal  model. 
Finally,  we  examine  the  variance  of  the  estimated  azimuth  for  various  array  geometries  in 
order  to  evaluate  the  proposed  triangular  local  arrays  with  various  baselines. 

The  above  results  enable  the  development  of  optimal  detectors  for  the  local  array  that  are 
based  on  frequency  domain  likelihood  arguments,  leading  to  statistics  that  are  monotone 
functions  of  the  usual  beam  power.  In  practice,  for  a  perfectly  correlated  signal  model 
the  optimal  test  statistic  is  the  ratio  of  the  beam  power  to  the  noise  power  and  has  the  F 
distribution.  In  the  decorrelated  signal  case,  which  may  be  common  when  sensors  in  the 
local  array  are  quite  far  apart,  we  obtain  the  likelihood  ratio  detector  as  a  weighted  linear 
combination  of  sample  cross-spectra. 

A  A 

Section  3  considers  optimal  methods  for  fusing  the  local  array  wave-number  vectors  9U . . . ,  0n 
into  an  overall  location  and  its  associated  uncertainty,  as  determined  by  90  and  95%  confi¬ 
dence  or  posterior  probability  ellipses.  We  consider  classical  nonlinear  least  squares  methods, 
with  ellipses  determined  assuming  the  variances  are  either  unknown  (Flinn,  1965),  unknown 
(Everndon,  1969)  or  unknown  with  a  specified  prior  distribution  (Jordan  and  Sverdrup, 
1981,  Bratt  and  Bache,  1988).  In  particular,  we  evaluate  the  performance  in  locating  the 
two  events  shown  in  Figure  2  as  a  function  of  several  array  geometries,  signal  to  noise  ratios 
and  detecting  configurations. 
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Section  2 

Local  Array  Detection  and  Estimation 


The  local  array  performance  can  be  evaluated  by  considering  the  model 

yi(t)  =  ei[t-7}(0)]  +  ni(t) 


(1) 


for  j  =  1, 2, . . . ,  N  sensors  recording  a  received  signal  Sj(t),  delayed  by  Tj(ff)  and  corrupted 
by  an  additive  noise  nj(t).  Note  that  9  —  is  a  vector  of  wave-number  parameters 

that  are  nonlinearly  related  to  velocity  and  azimuth  and  to  the  source  location  of  the  signal, 
say  x  =  (xi,xa)'.  If  the  coordinates  of  the  jth  sensor  are  denoted  by  rj  =  (r^i,r^)'  and  the 
signal  is  a  plane  wave  at  frequency  /,  measured  in  cycles  per  point,  we  may  take 


(2) 


as  giving  the  relation  between  the  time  delay  and  the  location  of  the  jth  sensor.  Furthermore, 


(3) 


and 

a=ta”~‘(t)  <4) 

give  the  relations  to  velocity  V  and  azimuth  a,  where  or  is  measured  clockwise  in  radians. 
Generally  we  convert  to  degrees  by  multiplying  by  180/tt  and  interpret  the  angle  obtained, 
beginning  with  0  degrees  north.  ||A||2  =  tr  AA*  will  denote  the  usual  squared  norm  of  the 
matrix  A,  where  *  denotes  the  complex  conjugate  transpose  of  the  matrix. 

For  assumptions,  we  will  always  take  the  noise  to  be  spatially  white  with  N  x  N  spectral 
matrix  P, ,( f)Ifr ,  where  In  denotes  the  N  x  N  identity  matrix.  The  signal  vector  st  = 
(«i(t), . . .  ,sw(<))'  is  assumed  to  be  stationary  with  unknown  NxN  spectral  matrix  E(/)  = 
{aii (/) >  *>  j  ~  1|  •  •  • » N}  in  the  general  decorrelated  signal  case,  with  squared  signal  coherence 


M/)  l2 

*«(/)*;;(/)’ 


which  might  possibly  form  the  basis  for  a  signal  model.  As  is  customary  in  these  kinds  of 
problems,  the  signal  and  noise  are  assumed  to  be  completely  uncorrelated. 


A  simplification  of  the  above  model  is  the  case  of  perfect  signal  correlation,  where  one  might 
as  well  assume  the  signal  to  be  identical  and  random  on  each  sensor,  i.e. 


Vj(t)  =  s[t  -  Tj(0))  +  nj(t), 


(5) 
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with  the  added  noise  being  uncorrelated  between  sensors  and  with  identical  spectrum  Pn(f) 
on  each  sensor.  As  will  be  seen  later,  this  simplifies  the  analysis  since  we  may  now  assume 
that  the  signal  spectrum  is  P,{f)  on  each  channel  and  that  the  signal  spectral  matrix  is  of 
the  form  P,(f)Iit . 

It  is  convenient  to  consider  first  the  maximum  likelihood  estimation  of  the  wave-number 
parameters  for  the  two  signal  models  (1)  and  (5)  under  the  assumption  that  the  signal  and 
noise  spectral  matrices  are  known.  Then,  we  move  to  estimating  the  signal  and  noise  spectral 
matrices,  conditionally  on  a  fixed  wave-number,  by  maximum  likelihood.  The  optimum 
detector  is  based  on  the  previous  two  results  and  the  distribution  theory  of  the  F-statistic 
obtained  allows  a  quantitative  assessment  for  the  significance  of  the  detected  signal. 


2.1  Estimation  of  Wave-Number  Parameters. 


The  estimation  of  the  wave-number  vector  0  is  most  easily  approached  using  the  Whit¬ 
tle  likelihood,  expressed  in  terms  of  the  probability  approximate  large  sample  probability 
distribution  of  the  discrete  Fourier  transform  (DFT)  of  the  observed  data,  namely 

=  (6) 

t=0 

where  ft  =  —  1,  ...,L  are  a  collection  of  frequencies  in  the  neighborhood  of  some 

center  frequency  /,  with  lT  a  sequence  of  adjacent  integers,  chosen  so  that  the  collected  IrfT 
form  a  band  that  encloses  the  desired  center  frequency.  The  Nxl  vectors  Y(  =  (Yu, . . . ,  Ystf 
are  assumed  to  have  a  common  multivariate  complex  normal  distribution  with  zero  mean 
and  complex  covariance  matrix 


£„(/;«)  =  P.(/M9)z(9)-+W)/„,  (7) 


where 


z(0)  =  (e7*ir,iB,...,e**ir'"0y 

denotes  the  Fourier  transform  of  the  vector  of  delays.  Details  involving  the  maximizing  the 
of  the  Whittle  likelihood  for  the  two  cases  below  can  be  found  in  Shumway  et  al  (1999). 


Case  i:  Perfectly  Correlated  Signals 


In  the  case  of  the  perfectly  correlated  signal  model  (5),  the  spectral  matrix  will  be  of  the 
form  (7)  and  the  log  likelihood  will  be  a  monotone  function  of  the  beam  power 


-2 w»r} 


(8) 


which  is  the  usual  observed  power  along  a  direction  corresponding  to  the  wave-number  0. 
Note  that  this  is  the  usual  beam  power  and  is  often  realized  by  filtering  the  received  data 
using  a  bandpass  filter  and  then  performing  the  beam  operation  in  the  time  domain. 
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The  main  problem  with  using  (8)  directly  as  a  signal  detector  is  that  its  performance  depends 
on  the  unknown  signal  and  noise  spectra,  P,(f)  and  Pn{f)  under  the  hypothesis  that  a  signal 
is  present  and  on  Pn(f)  when  the  signal  is  absent.  We  note  that  the  approximate  distribution 
of  the  test  statistic  is 

where  the  notation  denotes  a  chi-squared  random  variable  with  m  degrees  of  freedom. 
Hence,  the  detection  results  will  depend  upon  the  number  of  sensors  and  the  the  unknown 
signal  and  noise  parameters.  The  maximum  likelihood  estimators  for  these  nuisance  param¬ 
eters  are  given  in  Section  2.2  and  we  note  that  there  will  be  an  optimal  P-statistic  in  Section 
2.3  with  a  distribution  under  noise  alone  that  does  not  depend  on  the  unknown  signal  and 
noise  spectra. 


In  order  to  estimate  the  wave-number  vector,  it  is  conventional  to  plot  (8)  as  a  function  of 
0  and  to  choose  0  as  the  maximizer  of  B(6)  .  The  velocity  and  azimuth  then  follow  from 
(3)  and  (4).  The  large-sample  variances  and  covariances  of  the  maximizers  of  (8)  have  been 
derived  by  Shumway  et  al  (1999)  using  the  Cram6r-Rao  lower  bound.  In  this  simple  case, 
we  obtain  the  large-sample  covariance  matrices  as  a  function  of  the  covariance  matrix  of  the 
array  coordinates 

R=  jf!2(rj-f)(ri-rY  (9) 

iV  i= i 


and  the  signal-to-noise  ratio 


.  P.U) 
-Pn(f)' 


(10) 


A 

We  have  that  0  will  be  approximately  normal  with  mean  equal  to  the  true  mean  and  ap¬ 
proximate  covariance  matrix 


cov(0) 


1  1  1 
2(2tt)3  L  rN 


(11) 


Then,  defining  the  vectors  0  —  (0i,  62)'  and  9  =  (02,  —  0i)',  we  have  the  covariance  matrix  of 
velocity  and  azimuth  given  as 


cov 


ft# 


1  1  1 
2(2  n)*LrN 


(12) 


where 


(  ptfAO  -f\\0\\0'A9\ 


r(4)  = 


I  -S\\0\\«AO 


\\ofi/Ao  I 


(13) 


denotes  a  matrix  function  of  the  matrix  A  here  and  in  the  sequel.  This  exhibits  nicely 
the  dependence  of  the  variances  of  azimuth  and  velocity  on  the  geometry  of  the  array,  as 
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embodied  in  the  covariance  matrix  of  locations  R.  Note  the  variances  also  depend  on  the 
signal  to  noise  ratio  r,  so  that  we  still  need  an  estimator  for  this  quantity.  The  simplicity  of 
the  formulation  in  terms  of  the  wave-number  parameters  is  also  appealing  and  it  is  the  form 
of  (11)  as  opposed  to  (12)  that  suggests  the  use  of  wave-number  in  the  location  procedure.  In 
the  Section  2.3,  we  consider  estimating  the  signal  and  noise  spectra  using  a  spectral  analysis 
of  variance. 

Case  ii:  Signal  Decorrelation 

Introducing  signal  decorrelation  as  a  complication  in  the  model  also  complicates  the  estimar 
tion  of  the  wave-number  and  velocity  and  azimuth  parameters.  The  likelihood  corresponding 
to  the  general  signal  model  (7)  with  signal  decorrelation  leads  to  maximizing  a  generalized 
beam  power  of  the  form 

GB(6)  =  -^21 E  Yjt  Yu  (14) 

j,k  t 

where  the  weights  Cjk{f)  are  elements  of  the  matrix 

C(/)  =  [/k  +  P„(/)£-1(/)]-1,  (15) 

depending  specifically  on  the  spectral  matrix  of  the  signal.  We  discuss  a  procedure  for 
estimating  this  matrix  in  the  next  section.  Following  through  the  asymptotic  theory  for  the 
maximizer  0G  of  GD(0)  leads  to  a  variance  covariance  matrix  of  the  form 

°>v  fe)  «  (2*)/^C-'(/).  (16) 

where 

DU)  *  £ Cjk{f){rj  -  rOfo  -  rk)'akj(f)  (17) 

i* 

Again,  we  may  convert  into  the  covariance  matrix  for  velocity  and  azimuth,  obtaining 

COT(S)Ki45riPFr(D")’  (18) 

where  the  matrix  function  T(*)  has  been  defined  previously  in  (13). 

Because  of  the  increased  complexity  of  the  maximum  likelihood  estimator  for  wave-number 
in  the  decorrelated  case,  it  may  be  more  useful  to  evaluate  the  classical  beam  power  estimator 
for  this  case.  Shumway  et  al  (1999)  show  that  for  the  estimator  maximizing  the  beam  power, 
say  9,  it  follows  that 

cov  (0)  «  DilWDilt  (19) 
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where  Dx  is  just  D  in  (17),  with  Cj*  =  1  and 

7^2  =  E  {rj-nWt-rvfojflfpvkU) 

l27r'  ted'# 

+  2Pn(/)  £  fa  -  r*)fa  -  r*)W(/)  +  2P^(f)N^R.  (20) 

iJk# 

The  covariance  matrix  of  the  velocity  azimuth  pair  in  this  case  becomes 

“v(d)aj(4^rpFr(Dr‘H,Dr‘)-  (21) 

The  above  results  characterize  the  variances  and  covariances  of  the  usual  beam  power  under 
both  the  assumption  that  the  signals  are  perfectly  correlated  and  under  the  assumption 
that  there  is  a  coherence  that  declines,  possibly  as  a  function  of  distance.  The  next  section 
examines  these  quantities  as  a  function  of  array  geometry  at  a  fixed  signal  to  noise  ratio.  It 
should  be  noted  that  the  equations  given  in  this  section  provide  a  natural  set  of  estimated 

wave-numbers  from  each  of  n  possible  recording  arrays,  say  6i . On  and  a  set  of  estimated 

covariance  matrices,  say  cov  fa),. .  .,cov  fa,),  associated  with  these  estimated  parameters. 
The  results  enable  studying  the  problem  of  designing  arrays  as  in  the  next  section  as  well  as 
the  problem  of  fusing  single  array  parameters  into  an  optimal  location  as  in  Section  3. 

2.2  Theoretical  Performance  of  Single  Arrays. 

It  is  of  interest  to  determine  the  array  configuration  that  might  be  nearly  optimal  for  de¬ 
tection  and  azimuth  estimation  of  surface  and  underground  explosions.  We  note  that  the 
analysis  is  complicated  by  the  fact  that  the  signal  correlation  will  decrease  as  the  separation 
between  the  sensors  increases.  Blandford  (1996)  has  made  extrapolations  for  coherence  as  a 
function  of  distance,  using  a  model  of  Mack  and  Flinn  (1971),  for  the  distance  ranges  implied 
by  an  array  that  will  have  three  elements  arranged  in  a  triangle  with  vertices  separated  by  d 
kilometers  and  a  center  element.  Generally  speaking,  coherence  is  reasonably  high  at  low  fre¬ 
quencies  that  are  less  than,  say  .5  Hz.,  but  that  predicted  coherence  goes  down  significantly 
at  higher  frequencies,  particularly  for  the  separation  of  around  1  km  that  predominates  for 
the  triangular  array  under  consideration. 

It  is  interesting  to  compute  the  azimuthal  standard  deviations  implied  by  the  Mack-Flinn 
coherence  model  for  various  intersensor  separations  on  the  simple  triangular  array  with 
a  center  element.  The  approach  taken  here  assumes  that  either  the  values  maximizing 
the  beam  or  the  generalized  beam  will  be  used,  with  the  variances  computed  from  the 
appropriate  asymptotic  expressions  (18)  and  (21)  respectively.  Figure  3  shows  the  results  at 
various  distance  ranges  and  periods.  It  is  clear  that  the  beam  and  generalized  beam  behave 
similarly,  except  at  the  longer  distance  ranges.  The  signal  to  noise  ratio  r  was  taken  as 
.752,  which  is  regarded  as  sufficient  for  an  analyst  to  declare  a  detection.  Bandwidths  were 
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.01,  .02,  .04  and  .2  Hz  respectively  for  cases  (a)-(d)  in  Figure  3.  Azimuthal  uncertainties  are 
larger  for  longer  periods  and  actually  get  quite  small  for  the  higher  frequencies.  One  second 
period  standard  deviations  remain  below  1  degree  for  intersensor  distances  less  than  one  km. 
It  is  interesting  to  note  that  the  minimum  variance  is  obtained  for  separation  distances  that 
are  exactly  equal  to  the  period,  i.e.  20  km  for  a  20  sec  period,  10  km  for  a  10  second  period, 
5  km  for  a  5  second  period  and  1  km  for  a  1  second  period. 


Figure  3. 


(a)  20  sec.  Period  (b)  10  sec.  Period 


Estimated  large  sample  standard  deviation  of  azimuth  estimators  as  a  func¬ 
tion  of  baseline  distance  using  the  Mack-Flinn  theoretical  coherence  model  for 
decorrelation  with  distance  with  the  solid  lines  for  the  weighted  beamforming 
estimator  and  the  dotted  lines  denoting  simple  beamforming. 


Blandford  (1996)  also  increases  the  number  of  elements,  considering  an  N  =  7  element  array 
consisting  of  the  array  above  plus  a  small  inverted  interior  triangle,  with  .2  km  on  a  side. 
With  a  signal  to  noise  ratio  of  (.75 )2,  standard  deviations  were  in  the  .5-2.0  degree  range 
for  periods  less  than  5  seconds.  However,  there  is  little  improvement  for  10  and  20  second 
periods.  Blandford  (1996)  summarizes  the  computation  by  noting  that  for  the  5  second 
period  signal,  which  is  expected  from  a  1  kt  atmospheric  nuclear  explosion,  the  1  kilometer 
array  would  have  an  azimuth  estimation  error  approximately  equal  to  the  best  historically 
observed  residual  error.  The  conclusions  are  somewhat  preliminary  because  the  underlying 
data  is  available  only  for  much  long  periods  and  small  inter-sensor  distances.  Analysis  of 
additional  signal  data  at  more  appropriate  periods  and  spacings  is  of  critical  importance. 
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2.3  Estimation  of  Signal  and  Noise  Parameters. 

In  the  previous  sections,  we  have  regarded  the  signal  and  noise  spectra,  P,(/),  Pn(f)  in  the 
perfectly  correlated  case  to  be  known.  In  the  decorrelated  case,  we  need  to  know  the  entire 
N  x  N  spectral  signal  matrix  £(/).  It  turns  out  that  maximum  likelihood  estimators  are 
available  for  these  nuisance  parameters  and  these  estimators  are  summarized  in  the  following 
subsections. 


Case  i:  Perfectly  Correlated  Signals 


For  the  perfectly  correlated  case  the  log  likelihood  reduces  to 


In  L(0,  P„  Pn)  oc  -NL In  P»  -  L  ln(l  +  rN )  + 


m 

P„(AT  +  l/r) 


and  may  be  maximized  directly,  as  a  function  of  the  signal  to  noise  ratio  r  and  Pn. 
obtain,  for  a  fixed  value  of  0 , 


SSE(0 ) 
(JV-l)L’ 


We 

(23) 


where  SSE(0)  is  an  error  sum  of  squares  about  the  fitted  signal  and  may  be  written  as 


s5^(»)=f:E 

i=i  i=i  ;v 

The  signal  spectrum  is  estimated  by 


(24) 


P.(»)  =  ~  NLPM)  (25) 

We  may  then  apply  maximization  alternatively,  first  for  0  over  D(0)  and  then  solve  for 
P,(9)  and  Pn(0)  using  (23)  and  (25).  A  new  value  of  0  may  be  estimated  by  maximizing 
the  log  likelihood  (22)  and  the  obvious  iterative  pattern  will  converge  by  standard  alternate 
maximization  results  (see  Meng  and  Rubin,  1991).  Alternatively,  one  may  simply  stop  after 
the  first  step. 

Case  ii:  Signal  Decorrelation 


For  the  multivariate  signal  model  with  signal  spectral  matrix  E(/),  there  are  no  simple  es¬ 
timators  but  we  may  still  derive  a  procedure  for  maximizing  the  log  likelihood  using  either 
a  Newton-Raphson  approach  and  the  vector  scores  and  information  matrix  or  the  EM  al¬ 
gorithm  on  the  complete  data  log  likelihood.  Now,  regarding  the  generating  model  in  the 
frequency  domain  as 

Yji  —  Zj(0)Sji  +  Nji,  (26) 

we  may  write  the  vector  version  as 


Y(  =  Z{0)Si  +  Nt, 


(27) 


10 


where  Z{0)  =  diag(zi(0), . . . ,  zN(0))  is  a  diagonal  and  St  =  (Su>  5«, . . . ,  SW)'  is  the  random 
frequency  domain  signal  vector  with  spectral  matrix  £(/).  Considering  the  complete  data, 
in  the  sense  of  Dempster  et  al  (1977),  as  the  signal  vector  above  plus  the  noise  vectors 
Nu  •  • . ,  Nl,  we  obtain 

lnl(0,£,Pn)  *  -LlnlSl-E^S^ 

1=1 

-  IWlnP„-i-£||W,||2 3 4  (28) 

for  the  complete-data  log  likelihood.  To  maximize  the  original  incomplete  data  log  likelihood, 
we  consider  applying  the  Expected  Conditional  Maximization  (ECM)  algorithm  of  Meng 
and  Rubin  (1993).  In  this  version  of  the  EM  algorithm,  we  apply  the  expectation  and 
maximization  steps  for  a  given  0,  where  0  is  the  maximizer  for  a  fixed  £,  Pn.  Note  that  the 
EM  algorithm  of  Dempster  Laird  and  Rubin  (1977)  is  a  rather  convenient  way  of  handling 
maximum  likelihood  estimation  in  multidimensional  random  signal  context  (see  Shumway, 
1988).  Now  applying  the  E-step  of  the  EM  algorithm  to  the  complete  data  log  likelihood 
above  for  a  fixed  0  gives 

E  =  L-'£(s,  (9)5;  (9) +  £(«))  (29) 

for  the  next  iterate  of  the  signal  spectral  matrix,  where 


S,(9)  =  (/  +  PnS-l)',Z’(9)K/ 


GUiU  j 

E(8)  =  P„(/  +  P„ir1)’  , 

using  results  from  Shumway  et  al  (1999).  The  estimator  for  the  noise  spectrum  will  be 
P„  =  (NL)-'  £(r,  -  Z(9)S,(«)||2  +  tr(Z(8)£(9)Z‘(9))) 

To  apply  the  iterative  ECM  algorithm,  one  might  use  the  following  sequence  of  steps. 


1.  Determine  an  initial  estimator  for  0  as  the  maximizer  of  the  beam  power  B(0)  as  given 
by  (8). 

2.  Compute  initial  estimators  for  Pn  and  £  =  PJ,  assuming  no  signal  correlation,  using 
(23)  and  (25). 

3.  Update  5<(0)  and  £(0)  from  (30)  and  (31). 

4.  Update  £  and  Pn  using  (29)  and  (32). 
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5.  Repeat  steps  3.  and  4.  to  convergence 

6.  Update  0  by  finding  the  maximizer  of  GB(0)  in  (29). 

7.  Return  to  3. 


The  above  procedure  is  rather  involved  and  it  is  often  the  case  that  one  can  estimate  the 
signal  spectral  matrix  from  prior  data,  assuming  that  the  signal  spectral  matrix  in  known 
in  the  generalized  beamformed  estimator  GB(8).  This  is  illustrated  in  Section  2.2  where  we 
assigned  a  functional  form  to  Ip,-*),  where  Ojkif)  =  Pt(f)\pjk\  based  on  a  theoretical  relation 
between  coherence  and  distance. 


2.4  Signal  Detection. 

Detecting  a  signal  depends  on  developing  a  test  statistic  that  does  not  depend  on  any  nuisance 
parameters  under  noise  alone  and  is  a  convenient  function  of  the  signal  spectra  under  the 
signal  hypothesis.  We  indicate  below  how  to  use  the  results  of  the  previous  section  to  develop 
the  F-statistic  as  an  optimal  detector  in  the  case  of  a  perfectly  correlated  signal  and  give 
some  results  on  the  expected  performance  of  the  beam  under  the  decorrelated  signal  case. 

The  results  enable  a  prediction  of  the  performance  of  each  local  array  as  a  function  of  signal 
to  noise  ratio  in  the  perfectly  correlated  case  or  as  a  function  of  the  spectral  matrix  of  the 
signal  in  the  decorrelated  signal  case. 

Case  i:  Perfectly  Correlated  Signals 

Shumway  et  al  (1999)  establish  that  the  test  statistic 


(*-l) 

N 


B0) 

SSEifi) 


(33) 


has  asymptotically  a  central  F-distribution  with  with  2 L  and  2 L{N  -  1)  degrees  of  freedom 
when  there  is  no  signal,  where  B{0)  is  the  beam  power  defined  in  (8)  and  SSE(0)  is  the  error 
power  defined  in  (24).  The  F-statistic  is  evaluated  at  the  estimated  0  and  the  approximate 
distribution  is  given  by 

F(0)  ~  (1  +  NrfFziflurt-i)  (34) 

under  the  hypothesis  that  the  signal  is  present,  giving  a  computational  method  for  estab¬ 
lishing  the  false  alarm  rate  and  detection  probabilities  of  the  detector  as  a  function  of  the 
signal  to  noise  ratio  r  =  P,/Pn. 

Case  ii:  Signal  Decorrelation 

In  the  case  of  a  decorrelated  signal,  the  optimal  detector  becomes  the  likelihood  ratio  detec¬ 
tor,  which  will  be  a  monotone  function  of  the  generalized  beam  G(0 ).  One  can  follow  this 
argument  through  as  in  Shumway  et  al  (1999)  but  we  simply  mention  here  that  it  may  be 
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more  realistic  to  present  results  involving  the  simple  beamforming  statistic,  written  in  the 
form  L 

B(l)  =  '£PWt\1 

fc=l 

Now,  for  the  general  correlated  signal  model,  we  have  that  Yt  have  mean  zero  and  covariance 
matrix 

Ev(/)  =  Z(»)E(/)Z'(«)  +  />.(/)/„.  (35) 

It  follows  that  /  . 

B(«)  ~ 

which  simplifies,  on  substituting  on  the  right  hand  side  from  (35)  to 


B(«)  ~  +  NPM))ti»L, 


(36) 


and  enables  the  prediction  of  the  performance  of  the  beam  as  a  function  of  the  signal  spectral 
matrix  and  noise  spectra.  Note  that  we  have  replaced  the  wave-number  by  an  estimator  0 
converging  to  0  in  probability. 


2.5  Empirical  Results. 

We  look  first  at  a  simple  example  involving  the  three  sensors  in  Figure  1,  that  recorded  the 
event  Tanana  at  a  small  triangular  array  on  Palmyra  Island.  To  provide  a  baseline  example, 
we  took  2048  points  (about  200  seconds)  of  data  and  re-aligned  the  signal  to  correspond 
to  an  approximate  velocity  of  .3  km/sec  and  an  approximate  azimuth  of  225  degrees.  The 
array  geometry  was  a  triangle  with  1  km  sides  at  locations  rj  =  (0,  .577)',  r2  =  (.5,  —.289)' 
and  f3  as  (-.5,  -.289)',  leading  to  a  diagonal  covariance  matrix  R  in  (9)  that  has  1/6  on  the 
diagonals.  The  re-aligned  signals  at  sensors  A,B  and  C  are  shown  in  Figure  1;  the  signal 
spectrum  was  centered  at  .044  Hz  (cycles  per  second),  i.e.,  at  a  period  of  about  23  sec.  The 
coherence  between  sensors  was  reasonably  high,  in  this  case,  with  values  .86,  .83  and  .82  for 
the  three  inter-sensor  pairs.  Hence,  there  is  reason  to  believe  that  random  univariate  signal 
model  (5)  will  work  reasonably  well. 

In  order  for  the  large  L  asymptotics  to  have  a  chance,  we  first  chose  a  band  broad  enough  to 
include  the  entire  signal,  with  L=17  frequency  ordinates  in  the  band  running  from  .005  to 
.08  Hz.  The  estimated  signal-to-noise  ratio  Pt/Pn  —  1  /r,  in  this  band,  using  the  maximum 
likelihood  estimators  (23)  and  (25),  was  3.16.  For  the  broad  band,  the  estimator  for  the 
velocity  was  .26(.03)  km/sec,  using  (8)  and  (4)  to  get  the  estimator  and  (12)  to  get  the 
standard  deviation  (in  parentheses).  This  compares  to  the  velocity  of  .30  km/sec  that  we 
that  input  by  lining  up  the  largest  minima  of  the  signal  on  each  sensor  with  this  velocity.  The 
azimuth  estimator  is  225(7)  degrees,  which  is  right  on  the  known  azimuth.  The  standard 
deviation  implies  an  approximate  95%  confidence  interval  from  211  to  239  degrees.  For 
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(a)  Skatcft  of  F-Statiatlc< narrow  band) 


-05  -05 


(b)  Contour  of  F*Stad«tlca(narrow  band) 


Figure  4. 


(c)  Sketch  of  F-Statistic(broad  band)  (d)  Contour  of  F-Sutiitfcs (broad  band) 


Relief  (left  column)  and  contour  (right  column)  plots  of  F-statistics  for  detect¬ 
ing  infrasonic  signal  at  two  bandwidths;  narrow  band  F  shows  peak  at  227(10) 
degrees  and  velocity  c  =  .24(.04)  km/sec. 


comparison,  we  tried  a  narrow  band  estimator,  based  on  L=3  frequency  ordinates,  in  the 
band  running  from  approximately  .040  to  .048  Hz.  The  signal-to-noise  ratio  in  this  band  was 
7.36  and  the  estimated  velocity  of  .24(.05)  km/sec  is  still  slightly  off.  Again,  the  azimuth 
estimator  was  227(10)  degrees  with  a  slightly  larger  uncertainty.  The  increase  in  the  signal- 
to-noise  ratio  in  the  narrower  band  does  not  completely  compensate  for  the  loss  in  bandwidth. 

Figure  4  shows  the  relief  plots  (a)  and  (c)  and  contour  plots  (b)  and  (d)  of  the  F-detector  given 
by  the  F-statistic  (41).  The  detector  is  plotted  as  a  function  of  the  wave  number  coordinates 
—.5  <  0i,02  <  *5,  where  the  slowest  velocity  in  the  plot  corresponds  to  the  velocity  on  one 
of  the  edges,  e.g.  $i  =  =  *5,  or  .062  km/sec,  using  (4).  Very  high  velocities  occur  near 
the  center  0i  =  62  —  0,  which  corresponds  to  an  infinite  velocity.  The  azimuth  can  be  seen 
visually  as  the  angle  made  by  a  line  drawn  from  the  center  to  the  maximum;  here,  it  is  225 
degrees.  Note  that  the  asymptotics  for  this  will  be  less  sensitive  to  large  L  results,  although 
one  will  still  need  to  have  the  sample  length  T  large.  The  plotted  F-Statistic,  for  the  narrow 
band  case  has  2 L  =  2(3)  =  6  and  2 L[N  —  1)  =  6(2)  =  12  degrees  of  freedom.  Note  that 
■Fe,w(*01)  =  4.82  and  F#,i2(*001)  as  8.38  so  that  the  .01  and  .001  significance  points  are 
exceeded  by  all  velocities  and  azimuths  within  the  first  two  contours  of  the  top  plot.  The 
plotted  F-Statistic,  for  the  broad  band  case  has  2 L  =  2(17)  =  34  and  2L(N-1)  =  34(2)  =  68 
degrees  of  freedom.  Note  that  F34i68(.01)  =  1.95  and  FWiC8(.001)  =  2.42  so  that  the  .01  and 
.001  significance  points  are  exceeded  by  all  velocities  and  azimuths  within  the  first  two 
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Figure  5.  Relief  (left)  and  contour  (right)  plots  of  beam  power  and  maximum  likelihood 
correlated  beam  power  detectors.  Beam  power  shows  peak  at  225(7)  degrees 
and  velocity  V  =  .26(.03)  km/sec.,  correlated  detector  shows  peaks  at  225(6) 
degrees  and  .26(.03)  km/sec. 


contours  of  the  bottom  plot. 

It  is  also  of  interest  to  determine  the  estimated  velocities  and  azimuths  and  their  uncertainties 
under  the  coherence  structure  estimated  at  the  array.  This  requires  estimating  the  variance 
of  the  beamformed  estimator  maximizing  B(0)  under  the  assumption  that  there  is  less  than 
perfect  coherence  using  (21).  This  can  be  compared  to  the  variance  of  generalized  beam 
maximizing  GB(6),  given  by  (18).  Both  equations  require  estimating  the  N  x  N  signal 
spectral  matrix  £(/)  and  the  noise  spectrum  Pn(/).  In  the  case  of  the  0  maximizing  the 
beamforming  estimator,  we  used  (29)-(32)  and  the  EM  steps  1  to  5  to  obtain  the  maximum 
likelihood  estimators  for  E (/)  and  Pn(/),  with  0  fixed  at  0.  For  the  case  of  0D,  the  maximizer 
of  the  general  beam,  we  simply  used  the  ECM  algorithm,  following  Steps  1-7  as  given.  Figure 
5  shows  the  results  using  a  slightly  broader  band  (L=19)  and  we  note  that  the  estimators 
and  their  standard  errors  are  nearly  the  same,  a  reflection,  no  doubt  of  the  relatively  high 
empirical  coherence  between  the  three  sensors  in  this  case. 
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Section  3 

Global  Location  Capabilities 


Integrating  or  fusing  data  from  the  single-array  sources  into  a  best  overall  location,  with 
an  uncertainty  region  provided,  will  be  an  important  aspect  of  evaluating  the  predicted 
performance  of  the  IMS  network  shown  in  Figure  2.  Hence,  we  consider  a  methodology  for 
using  the  information  developed  in  the  detection  and  estimation  portions  of  this  proposal 
for  estimating  the  location  vector  x  =  (*j,  x2)'-  In  general,  we  propose  a  model  of  the  form 


8k  =  Bk(x)  +  ek,  (37) 

where  d*  =  (^i*, A:  =  are  the  estimated  wave-number  vectors,  as  computed 

from  maximizing  the  beam  power  (8),  or  F-statistic  (33)  at  the  kth  array  and 


hi*)  = 


fk  X  —  Ck 


gives  the  theoretical  connection  between  the  wave-number  parameters  and  location.  In  (38), 
fk  is  the  center  frequency,  Vk  is  velocity  and  ck  -  (ci*,  c^)'  denotes  the  coordinates  of  the  Jfctt 
array.  It  can  generally  be  assumed  that  velocity  is  known  or  can  be  inferred  from  the  wave- 
number  plot.  It  should  be  noted  that  there  are  often  separate  phases  at  the  same  array  that 
may  have  their  own  estimated  Ok  and  these  are  included  under  those  that  possibly  contribute 
to  the  model  (38).  Under  certain  conditions,  there  may  also  be  estimators  of  travel  times 
tk,  where  f*(x)  —  ||x  —  ck\\/Vk,  perhaps  from  cross-correlation  or  other  means,  that  could  be 
added  to  the  observations  on  location  x.  This  general  framework  also  does  not  preclude  the 
possibility  of  adding  seismic  travel  times  or  hyroacoustic  information  to  the  stack. 

A  possible  assumption  for  the  bivanate  error  terms  in  (37)  is  that  they  are  independent  and 
identically  distributed  with  mean  zero  and  2x2  covariance  matrix 


cov  cjt  =  <r2£(0*),  (39) 

where  the  matrix  £($*)  comes  from  (11),  specialized  to  the  kth  array.  Note  that  the  com¬ 
ponents  of  (11)  will  vary  according  to  the  array  size  and  geometry,  signal-to-noise  ratio  and 
the  time-bandwidth  product.  The  scaling  variance  <r2  is  to  account  for  additional  variability 
from  geophysical  sources  or  from  the  observed  error  in  a  particular  event  location.  We  have 
used  the  scaling  rather  an  additive  model  for  the  geophysical  error  expansion  to  simplify  the 
Bayesian  computation.  If  there  are  consistent  biases  associated  with  particular  regions  or 
subsets  of  arrays,  constant  correction  terms  can  be  added  to  the  defining  Equation  (37).  If  a 
number  of  events  are  available,  the  correction  terms  may  even  be  estimated  by  least  squares 
using  consistent  source-receiver  pairs. 
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As  an  example  of  a  small  demonstration  set,  we  consider  Figure  2,  which  shows  14  IMS  array 
centers  c*  proposed  for  the  southwest  quadrant  of  the  world-wide  network.  For  purposes 
of  illustration,  we  also  show  two  hypothetical  events  that  might  generate  observed  wave- 
number  parameters  $*,  k  =  1, . . . ,  n  for  incorporation  into  the  location  model  defined  in  (37) 
and  (38).  In  Section  3.1,  we  discuss  an  approach  to  combining  or  fusing  the  wave-number 
parameters  into  an  overall  location,  based  on  a  given  fixed  set  of  recording  arrays  and  their 
characteristics,  as  defined  by  the  signal-to-noise  ratio  r,  the  time-bandwidth  product  BT 
and  the  array  geometry,  as  it  is  expressed  in  terms  of  R  and  N  in  (11). 

3.1  Fusion  Locations  Via  Classical  and  Bayesian  Methods. 

We  extend  the  classical  methods  first  to  the  case  where  we  observe  wave-number  parameters 
and  their  covariance  matrix  from  n  arrays  and  wish  to  combine  or  fuse  the  information  into 
an  overall  location.  The  nonlinear  model  (37)  and  (38)  can  be  treated  in  the  usual  way. 
That  is,  expand  0*(x)  around  some  initial  value,  say  x  =  x<>  and  write  a  linearization  as 

Bk  -  Bk  (x0)  =  Ak (x0)(x  -  x0)  +  ek,  (40) 


where 


■A»(z) = 


aok(x) 

dx 


(41) 


is  the  usual  2x2  matrix  of  partial  derivatives  of  0*(x).  Then,  stacking  the  n,  2  x  1  wave- 
number  vectors  and  minimizing  the  weighted  sum  of  squared  errors  can  be  done  by  succes¬ 
sively  estimating  P  —  x  —  xq.  This  leads  to 


x  =  xo  +  C  1  (xo)  53  Ak(xo)'  E  1  (0*)[0*  —  0*(xo)], 

jfe=i 


(42) 


where  n 

C(x„)  =  £At(*o)'£-,(«»M»(*o)-  (43) 

k= 1 

It  follows  that  the  estimated  covariance  matrix  of  the  final  estimator  is 

cov  x  =  cr2C~l(x).  (44) 

Equations  (42)  and  (43)  exhibit  the  fusion  estimators  at  each  stage  as  pooled  estimators  over 
the  »  arrays  as  long  as  the  variances  are  known.  We  may  also  develop  a  confidence  ellipse 
for  the  fusion  estimators  under  assumptions  (A),  (B)  and  (C),  as  given  below. 

A.  Variance  Known 

We  may  assume  that  the  variance  a2  is  known,  either  from  the  statistical  variances  of  the 
computed  wave-number  estimators  or  from  a  combination  of  factors  including  the  statistical 
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wave-number  variances.  In  this  case,  the  generalization  of  the  usual  chi-squared  ellipse 
considered  by  Evemden  (1969)  can  be  computed  from  the  fact  that 

(z  -  x)'C(x)(x -x)~  (45) 

where  ~  denotes  is  distributed  as  and  x2  denotes  a  chi-squared  distribution  with  2  degrees 
of  freedom.  Note  that  the  statistical  uncertainty  of  the  wave-number  estimators  is  already 
in  the  matrix  E(0*)  so  that  a  plausible  estimator  for  a7  in  the  absence  of  other  factors  might 
be  unity. 

B.  Variance  Unknown 

If  variances  are  known  only  up  to  the  constant  a7,  this  scaling  variance  may  be  estimated 
from  the  set  of  arrays  that  record  the  event.  For  the  Gaussian  case,  the  maximum  likelihood 
estimator  is  proportional  to  the  unbiased  estimator 

£(»*  -**(*))'  £"&)(»*  -*.(*))•  (46) 

This  case,  originally  considered  in  Flinn(1965),  leads  to  a  confidence  interval  based  on  the 
F-distribution,  namely 

(x  -  x)'C(x)(x  -x)  ~  2s2F2i2(n_1),  (47) 

where  F2,2(n-i)  denotes  the  F-distribution  with  2  and  2(n  -  1)  degrees  of  freedom. 


Figure  6.  Possible  prior  distributions  for  standard  deviations  of  measured  wave-number 
estimates. 


C.  Variance  Subject  to  Prior  Distribution 

It  is  often  the  case  that  it  is  unrealistic  to  assume  that  the  variance  is  known  exactly  because 
the  ellipse  defined  by  (45)  becomes  too  small.  For  a  small  number  of  arrays,  the  ellipse  based 
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on  the  F-statistic  (47)  is  often  much  too  large.  A  useful  compromise,  introduced  by  Jordan 
and  Sverdrup  (1981)and  continued  by  Bratt  and  Bache  (1988),  is  to  quantify  the  initial 
uncertainty  about  a1  by  assigning  it  a  prior  distribution  with  density  function  7r(cr2).  It 
is  convenient  to  use  the  inverted  chi-squared  distribution  with  parameters  m,  representing 
the  equivalent  sample  size  embodied  in  the  prior  information  and  0%,  representing  a  prior 
centering  value  for  the  variance.  For  the  form  of  the  density  function,  see  Anderson  (1984). 
Figure  6  plots  the  density  function  for  the  standard  deviation  a  for  a0  -  1  and  m  -  10, 30. 
We  note  that  the  two  values  put  the  standard  deviation  between  .4  and  2  for  m  =  10  and 
between  .6  and  1.6  for  m  =  30  For  a  fully  Bayesian  approach,  we  assume  a  non-informative 
prior  on  (-00  <  xi,x2  <  00)  for  the  location  x  and  compute  the  posterior  distribution,  given 
the  wave-number  observations,  as  a  bivariate  t-distribution  with  2  and  2(n  —  1)  +  m  degrees 
of  freedom.  The  posterior  estimator  for  the  variance  is 


n  _  2(n  -  l)s2  +  mol 
2(n  -  1)  +  m  ' 


(48) 


implying  that  the  best  approach  is  simply  to  pool  the  initial  variance  a2  and  the  sample 
variance  s2,  weighted  by  their  degrees  of  freedom.  The  quadratic  form  involving  the  location 
vector  x  in  the  multivariate  t  has  an  F-distribution,  making  the  95%  posterior  probability 
ellipse  for  the  location  expressible  as 


(x  -  x)'C(x)(x  -x)~  2<7'2F2t2(„_i)+m  (49) 


It  is  interesting  that  the  form  of  the  posterior  probability  ellipse  (49)  is  similar  to  (47)  but 
will  be  tighter  because  of  the  additional  degrees  of  freedom  for  the  F-statistic.  Hence,  the 
Bayesian  solution  represents  a  compromise  between  (45)  and  (47),  the  methods  of  (A)  and 
(B). 


3.2  Locations  Using  the  IMS  Array. 

Given  the  theoretical  developments  of  the  previous  sections,  it  is  now  possible  to  begin  to 
assess  the  potential  location  capabilities  of  the  reduced  IMS  network  shown  in  Figure  2.  It  is 
likely  that  there  will  be  relatively  few  stations  recording  each  event  so  we  will  make  compu¬ 
tations  based  on  the  Bayesian  approach  described  in  (C)  of  the  previous  section.  The  natural 
measure  of  error  in  the  estimated  location  is  the  confidence  ellipse  in  the  classical  methods 
(A)  and  (B)  and  the  posterior  probability  ellipse  for  the  Bayesian  approach  described  in  (C). 

For  the  small-scale  example  considered  here,  the  following  inputs  are  needed  to  estimate 
location  capability.  We  confine  this  particular  example  to  the  perfectly  correlated  case 
although  it  is  obvious  that  this  would  not  be  constraint  because  of  the  separate  analysis 
available  for  decorrelated  arrays. 
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Signal  to  Noise  Ratio 

The  signal  to  noise  ratio  r,  defined  in  (10)  as  the  ratio  of  the  signal  spectrum  to  the  noise  spec¬ 
trum,  is  critical  for  determining  the  covariance  input  (11)  which  is  needed  for  the  covariance 
matrix  E ($*)  in  the  posterior  probability  intervals.  The  signal  to  noise  ratio  also  determines 
the  detection  probability  from  (34),  when  the  signal  is  present;  this  determines  the  probable 
configuration  of  recording  stations.  The  signal  to  noise  ratio  can  be  estimated  using  (23) 
and  (25)  if  data  are  available  or  it  can  be  assumed  from  known  operating  conditions.  For 
the  example  given  here,  it  is  assumed  to  be  r  =  2, 3. 

Array  Geometry 

The  geometry  of  each  local  array  affects  the  detection  probability  (34)  and  the  estimated 
covariance  of  the  input  wave-number  parameter  E(0*)  in  (11),  essentially  through  the  number 
of  sensors,  JV,  and  the  covariance  matrix  of  the  sensor  locations  R.  In  the  example  given 
next,  we  took  two  array  geometries  (see  Blandford,  1996),  namely  a  triangle  with  a  1  km 
baseline  and  a  center  element  (N=4)  and  and  extended  version  that  adds  an  interior  inverted 
inner  triangle  with  baseline  distance  .2  km  (N=7). 

Time-Bandwidth  Product,  Frequency,  Velocity 

The  time-bandwidth  product  L  =  BT  affects  both  the  covariance  matrix  E(0*)  in  (11) 
and  the  detection  probability  (34).  Again  this  will  change  both  the  input  covariance  for  the 
location  computation  and  will  be  a  determinant  of  the  configuration  of  detecting  local  arrays. 
We  assumed  here  a  center  frequency  of  1  Hz  and  a  time-bandwidth  product  of  BT  =  17, 
which  represented  a  reasonable  compromise  from  our  earlier  work  (Shumway  et  all,  1999) 
with  a  Pacific  Islands  event.  Velocity  was  assumed  to  be  .3  km/sec  in  the  example  but  it 
could  be  estimated  for  a  given  event  at  a  given  station  using  (3). 

Prior  Distribution  of  Scaling  Variance 

As  an  approximate  prior  distribution  for  the  scaling  variance  a2,  we  can  examine  plots  like 
Figure  6  for  various  degrees  of  freedom,  m,  which  can  be  regarded  as  the  equivalent  sample 
size  that  would  have  been  required  to  produce  the  required  degree  of  accuracy  for  a1.  For 
purposes  of  illustration,  assume  that  we  expect  the  variance  to  be  =  1  with  some  spread 
about  the  value,  corresponding  tom  =  10.  This  is  probably  optimistic  as  one  might  expect 
scaling  due  to  geophysical  causes  to  be  somewhat  greater  than  unity  but  the  value  1  provides 
a  reasonable  illustrative  constant.  Figure  6  shows  the  range  of  values  that  could  be  expected 
with  this  choice  and  we  assumed  that  the  sample  variance  s2  =  1  would  have  been  computed 
using  (46). 

Using  the  above  specifications  gives  the  posterior  estimators  for  the  variance  from  (48)  with 
n  =  4  and  n  =  7  corresponding  to  the  two  local  array  configurations.  Assuming  the  true 
location  at  x  =  xq,  corresponding  to  the  two  known  locations  marked  in  Figure  2  gives 
reasonable  scenarios  for  event  location.  The  posterior  probability  ellipses  can  be  computed 
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Figure  7.  Posterior  probability  (90  and  95%)  ellipses  for  various  array  detection  config¬ 
urations  assuming  7-element  arrays,  perfect  signal  correlation  and  S/N  —  2. 


using  (49)  as  a  function  of  x,  evaluated  at  x  —  2q,  with  C(x o)  computed  from  (43). 

Figure  7  shows  the  resulting  posterior  probability  ellipses  for  the  extended  array  (n  =  7), 
assuming  some  typical  recording  configurations.  For  Event  1,  we  assumed  two  configurations 
of  the  closest  three  stations  (5,8,10)  detected  the  event.  Note  that  the  area  of  the  90%  ellipse 
decreased  by  about  20%  in  going  from  two  to  three  station  detection.  For  Event  2,  comparing 
n  as  4  station  detection  against  n  =  7  station  detection  reduced  the  area  of  the  ellipse  by 
about  30%. 

Table  1  gives  a  more  detailed  accounting  as  to  what  happens  for  various  scenarios  over  the 
recording  configurations  in  Figure  7.  Here,  we  see  30-40%  reductions  achieved  by  increasing 
the  signal  to  noise  ratio  from  2  to  3.  Reductions  achieved  by  adding  the  three  center  elements 
are  more  modest,  on  the  order  of  7-10%. 
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Table  1:  Areas  (km2)  of  90%  Posterior  Probability  Ellipses  for  Simple  and  Extended 

Triangular  Arrays 


Triangular 

Array 

Extended 

Array 

Stations/Event 

S/N  =  2 

CO 

II 

S/N  =  2 

S/N  =  3 

5,8/1 

5554 

3566 

5087 

3316 

4448 

2856 

4074 

2656 

4,12,13/2 

3129 

2002 

2858 

1863 

1,2,4,11,12,13/2 

2176 

1397 

1992 

1299 
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Section  4 
Discussion 


In  this  project  we  have  developed  an  integrated  approach  to  estimating  velocities  and  az¬ 
imuths  from  a  collection  of  local  arrays  and  then  fusing  the  data  into  Bayesian  locations  and 
their  associated  uncertainty  ellipses.  A  small-array  theory  is  developed  that  characterizes  the 
performance  of  local  optimal  detectors  under  signal  correlation  and  decorrelation  scenarios. 
We  compare  the  performance  of  maximum  likelihood  estimators  such  as  the  beam  power 
and  the  generalized  beam  power  as  a  function  of  array  geometry  and  signal  to  noise  ratio. 
Optimal  local-array  geometries  axe  suggested  that  are  relevant  to  the  problem  of  designing 
an  optimal  infrasound  array. 

Wave-number  estimators  along  with  estimated  variance  covariance  matrices  are  used  as  input 
to  study  the  size  and  orientation  of  90%  posterior  probability  ellipses  for  various  likely 
subsets  of  detecting  stations  within  the  global  infrasound  array  proposed  for  the  Prototype 
International  Data  Center  (IMS).  Adding  detecting  stations  decreased  the  size  of  the  90% 
ellipse  by  about  10-20%  per  added  station,  whereas  increasing  the  signal  to  noise  ratio  from  2 
to  3  decreased  the  size  of  the  ellipse  by  30-40%.  Adding  an  inner  triangle  to  the  conventional 
1  km  triangular  array  gave  more  modest  reductions  of  7-10%. 

Since  the  above  conclusions  are  tentative  and  only  apply  to  limited  array  configurations, 
any  recommendations  based  on  them  should  be  considered  as  preliminary.  The  solution 
for  a  fixed  configuration  of  recording  arrays  is  not  the  only  factor  of  interest  since  there 
will  be  multiple  possible  recording  configurations  possible  for  any  given  event.  That  is,  the 
capability  of  the  network  will  be  an  expected  value,  accumulated  over  possible  recording 
configurations,  weighted  by  the  probability  of  each  particular  configuration.  The  resulting 
network  capability  can  be  contoured  by  setting  an  event  at  each  location  on  the  world-wide 
grid  and  then  computing  the  probability  that  each  station  detects  the  event .  We  may  then 
compute  a  weighted  average  of  some  parameter  reflecting  the  location  capability.  For  this 
discussion,  we  assume  that  the  area  of  the  90%  ellipse  is  of  interest,  where  the  ellipse  may 
be  a  confidence  ellipse  under  classical  assumptions  or  a  posterior  probability  ellipse  under 
the  Bayesian  paradigm. 

This  might  proceed  in  a  manner  that  essentially  parallels  the  Networth  calculations  made 
by  Wirth  et  al  (1976)  in  the  case  of  seismic  arrays.  It  would  be  most  useful  to  develop  an 
expression  for  the  average  area  expected  over  all  possible  configurations  of  detecting  arrays. 
That  is,  define  a  detection  indicator  £>*  that  is  one  if  the  array  detects  and  is  zero  otherwise, 
for  the  full  set  of  possible  detecting  arrays,  say,  for  k  =  1,2,..., K  of  them.  If  we  define 
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(50) 


Pt  =  Pr{Dt  =  1},  the  joint  density  of  the  random  variables  DUD2 . DK  would  be 

p(d . . 

k=l 

The  probabilities  can  be  arbitrary  for  each  station,  based  on  recorded  data  or  on  the  esti¬ 
mated  probabilities  computed  from  the  result  (4).  Now,  for  any  given  event,  we  can  observe 
2*  possible  configurations  of  detecting  stations  and  each  one  of  them  will  give  a  predicted 
area.  We  can  multiply  each  configuration  probability  by  its  area  and  add  them  up  to  get 
the  average  areas  of  location  for  that  particular  source.  Alternately,  and  perhaps,  easier 
would  be  to  simulate  values  of  . . . ,  D*  repeatedly  and  simply  average  the  resulting  ar¬ 
eas.  In  either  case,  it  is  clear  that  a  relatively  small  computing  effort  will  yield  an  average 
predicted  90%  uncertainty  area  for  each  location.  Plotting  plotting  the  results  as  a  grid  of 
location  contours  on  the  map  would  give  a  more  detailed  index  of  potential  global  network 
performance. 
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ALEXANDRIA,  VA  22310-3398 
ATTN:  CTR,  L.  KLUCHKO 
ATTN:  TD,  DR  D.  LINGER 
ATTN:  TDANP.TRC 
ATTN:  TDC,  M.  SHORE 
ATTN:  TDCN,  R.  JIH 

HEADQUARTERS 

DEFENSE  THREAT  REDUCTION  AGENCY 
8725  JOHN  J.  KINGMAN  ROAD,  MS-6201 
FT  BELVOIR,  VA  22060-6201 
ATTN:  TDCS,  DR  J.  FOX 

DEFENSE  THREAT  REDUCTION  AGENCY 
ALBUQUERQUE  OPERATIONS 
1680  TEXAS  STREET,  SE 
KIRTLAND  AFB,  NM  87117-5669 
ATTN:  TDT-D,  DR  G.  BALADI 
ATTN:  TDTP,  DR  B.  RISTVET 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 

ACIS 
6Ti1  NHB 

WASHINGTON,  DC  20505 
ATTN:  DR  J.  FILSON 
ATTN:  L.  TURNBULL 
ATTN:  T.  MURPHY 

BATTELLE  MEMORIAL  INSTITUTE 
MUNITIONS  &  ORDNANCE  CTR 
505  KING  AVENUE 
COLUMBUS,  OH  43201-2693 
ATTN:  TACTEC 


BBN  CORPORATION 
1300  N.  17TH  STREET,  SUITE  1200 
ARLINGTON,  VA  22209 
ATTN:  H.  FARRELL 
ATTN:  J.PULLI 


BDM  CORPORATION  OF  SAUDI  ARABIA 
12150  EAST  MONUMENT  DRIVE.  SUITE  501 
FAIRFAX,  VA  22033-4050 

ATTN:  J.  STOCKTON,  JB  4C22 

CALIFORNIA,  UNIVERSITY  AT  SAN  DIEGO 
SCRIPPS  INSTITUTION  OF  OCEANOGRAPHY 
P.  O.  BOX  6049 
SAN  DIEGO,  CA  92166-6049 
ATTN:  C.  DEGROOT-HEDUN 
ATTN:  DR  M.  A.  H.  HEDLIN 
ATTN:  PROF  F.  VERNON 
ATTN:  PROF  J.  A.  ORCUTT 
ATTN:  PROF  J.  BERGER 

CENTER  FOR  MONITORING  RESEARCH 
1300  N,  17TH  STREET,  SUITE  1450 
ARLINGTON,  VA  22209 

ATTN:  DR  K.  L.  MCLAUGHLIN 
ATTN:  DR  R.  NORTH 
ATTN:  LIBRARIAN 
ATTN:  R.  A.  GUSTAFSON 
ATTN:  V.  RYABOY 

CTB  TREATY  MANAGER 
ROSSLYN  GATEWAY 
1901  N.  MOORE  STREET,  SUITE  609 
ARLINGTON,  VA  22209 

ATTN:  DR  R.  W.  ALEWINE,  III 

DGI 

3612  SURFWOOD  RD 
MALIBU,  CA  90265 

ATTN:  R,  F.HERBST 

ENSCO,  INC 
445  PINEDA  COURT 
MELBOURNE,  FL  3940 
ATTN:  DR  D.  TAYLOR 


ENSCO,  INC. 

P.  O.  BOX  1346 

SPRINGFIELD,  VA  22151-0346 
ATTN:  D.  BAUMGARDT 
ATTN:  Z.DER 
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GEOPHEX,  LTD. 

WESTON  GEOPHYSICAL 
325  WEST  MAIN  STREET 
NORTHBOROUGH,  MA  01532 
ATTN:  DR  D.  REITER 
ATTN:  MR  J.  LEWKOWIC2 

ITT  INDUSTRIES 
ITT  SYSTEMS  CORPORATION 
ATTN:  AODTRA/DTRIAC 
1680  TEXAS  STREET.SE 
KIRTLAND  AFB,  NM  87117-5669 
ATTN:  DTRIAC 
ATTN:  DTRIAC/DARE 


JAYCOR 

1410  SPRING  HILL  ROAD,  SUITE  300 
MCLEAN,  VA  22102 

ATTN:  DR  C.  P.  KNOWLES 
ATTN:  T.  J.  HANNIGAN 

JAYCOR,  INC. 

P.O.BOX  85154 
SAN  DIEGO,  CA  92138-5154 
ATTN:  M.TREADAWAY 
ATTN:  R.  STAHL 

LACHEL  &  ASSOCIATES  INC. 

P.  O.  BOX  5266 
GOLDEN  CO  80401 
ATTN:  D.  LACHEL 

LOGICON  INC. 

LOGICON  ADVANCED  TECHNOLOGY 
2100  WASHINGTON  BLVD 
ARLINGTON,  VA  22204-5704 
ATTN:  F.  DICKERSON 

LOGICON  R  AND  D  ASSOCIATES 
P.  O.  BOX  100 
PAHRUMP,  NV  89041 
ATTN:  J.LACOMB 

MAXWELL  LABORATORIES,  INC. 

S-CUBED  WASHINGTON  RESEARCH  OFFICE 
11800  SUNRISE  VALLEY  DRIVE,  SUITE  1212 
RESTON,  VA  22091 

ATTN:  DR.  T.  J.  BENNETT 
ATTN:  J.R.  MURPHY 


MAXWELL  TECHNOLOGIES  INC. 
SYSTEMS  DIVISION 
9210  SKY  PARK  COURT 
SAN  DIEGO,  CA  92123-4302 
ATTN:  DR  E.  PETERSON 
ATTN:  DR  J.  L.  STEVENS 
ATTN:  DR  G.E.  BAKER 


MISSION  RESEARCH  CORP. 

8560  CINDERBED  ROAD,  SUITE  700 
NEWINGTON,  VA  22122 
ATTN:  DR.  M.  FISK 


MULTIMAX,  INC 
1441  MC  CORMICK  DRIVE 
LANDOVER,  MD  20785 
ATTN:  DR  I.N.  GUPTA 
ATTN:  DRW. CHAN 
ATTN:  MS  L.  GRANT 


PACIFIC  NORTHWEST  LABORATORIES 
A  DIV  OF  BATTELLE  MEMORIAL  INST. 

P.  O.  BOX  999 
RICHLAND,  MA  99352 

ATTN:  TECHNICAL  STAFF,  MS  K5-12 

PACIFIC-SIERRA  RESEARCH  CORP. 
WASHINGTON  OPERATIONS 
1400  KEY  BOULEVARD,  SUITE  700 
ARLINGTON,  VA  22209 
ATTN:  N.L.  DUNCAN 


S-CUBED 

A  DIVISION  OF  MAXWELL  LABS,  INC. 

11800  SUNRISE  VALLEY  DRIVE,  SUITE  1212 
RESTON,  VA  22091 
ATTN:  J.  MURPHY 


SCIENCE  &  ENGINEERING  ASSOCIATES,  INC. 
7918  JONES  BRANCH  DRIVE,  SUITE  500 
MCLEAN,  VA  22102 
ATTN:  R.  BEATY 


SCIENCE  APPLICATION  INT'L  CORP. 
3309  NW  GOLDEN  PLACE 
SEATTLE,  WA  98117 
ATTN:  A.  RATLETON 


SCIENCE  APPLICATIONS  INTL  CORP. 
10260  CAMPUS  POINT  DRIVE 
SAN  DIEGO,  CA  92121-1578 
ATTN:  DR  T.  C.  BACHE,  JR. 
ATTN:  DR  T.  J.  SERENO,  JR. 
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SCIENCE  APPLICATIONS  INTL  CORP 
21 1 1  EISENHOWER  AVENUE,  SUITE  205 
ALEXANDRIA,  VA  22314 

ATTN:  DR  D.  PIEPENBURG 

SCIENCE  APPLICATIONS  INTL  CORP. 
GEOSPATIAL  DATA  DEVELOPMENT  DIVISION 
700  SOUTH  BABCOCK  STREET,  SUITE  300 
MELBOURNE,  FL  32901 
ATTN:  R.  BJURSTROM 

SCIENCE  APPLICATIONS  INTL  CORP. 

P.O.  BOX  1303 
MCLEAN,  VA  22102 
ATTN:  D.  BACON 

SOUTHERN  METHODIST  UNIVERSITY 
DEPT  OF  GEOLOGICAL  SCIENCE 
P.O.  BOX 395 
DALLAS,  TX  75275*0395 
ATTN:  DR  B.  STUMP 
ATTN:  E.  HERRIN 

SRI  INTERNATIONAL 
333  RAVENSWOOD  AVENUE 
MENLO  PARK,  CA  94025-3434 
ATTN:  D.  CURRAN 

ST  LOUIS  UNIVERSITY 
P.  O.  BOX  8148 
PIERRE  LACLEDE  STATION 
ST  LOUIS,  MO  63156*8148 
ATTN:  PROF  B.  HERRMANN 
ATTN:  PROF  B.  J.  MITCHELL 

TASC,  INC 

1101  WILSON  BOULEVARD,  SUITE  1500 
ARLINGTON,  VA  22209*2248 
ATTN:  J.A.MOSORA 

TEXAS,  UNIVERSITY  AT  AUSTIN 
P.O.  BOX 7726 
AUSTIN,  TX  78712 

ATTN:  C.  A.  FROELICH 


WOODWARD-CLYDE  CONSULTANTS 
566  EL  DORADO  STREET 
PASADENA,  CA  91109-3245 
ATTN:  DR  B.B.  WOODS 
ATTN:  DRC.K.  SAIKIA 

DEPARTMENT  OF  ENERGY 


BECHTEL  NEVADA 
P.  O.  BOX  809 

LOS  ALAMOS,  NM  87544-0809 
ATTN:  D.EILERS 


BECHTEL  NEVADA,  INC 

P.  O.  SOX  98521 

LAS  VEGAS,  NV  89193-8521 

ATTN:  D.  BARKER.  M/S  NLV053 

DEPARTMENT  OF  ENERGY 
1000  INDEPENDENCE  AVENUE  SW 
WASHINGTON.  DC  20585 

ATTN:  D.  WATKINS,  NN-42/GA007 
ATTN:  E.  MANAK,  NN-30/FORS 
ATTN:  E.  STOVER 
ATTN:  S.  RUDNICK/NN-20 

UNIVERSITY  OF  CALIFORNIA 
LAWRENCE  LIVERMORE  NATIONAL  LAB 
P.  O.  BOX  808 

LIVERMORE,  CA  94551*9900 
ATTN:  K.  NAKANISHI 
ATTN:  W.  J.HANNON,  JR,  MS,  L-103 
ATTN:  F.  HEUZE.  MS,  L-200 
ATTN:  L.  GLENN.  MS,  L-200 
ATTN:  DR  J.  2UCCA,  MS,  L-205 
ATTN:  M.  DENNY,  MS,  L-205 
ATTN:  TECHNICAL  STAFF,  MS,  L-200 
ATTN:  TECHNICAL  STAFF,  MS,  L-208 
ATTN:  TECHNICAL  STAFF,  MS,  L-205 

LOS  ALAMOS  NATIONAL  LABORATORY 

P.O.  BOX  1663 

MAIL  STOP  G-733 

LOS  ALAMOS,  NM  87545 

ATTN:  F.  CHAVEZ,  MS-D460 
ATTN:  D.  WESTERVELT,  MS-A112 
ATTN:  D.  STEEDMAN,  MS-F607 
ATTN:  TECHNICAL  STAFF,  MS-C335 
ATTN:  TECHNICAL  STAFF.  MS-D460 
ATTN:  TECHNICAL  STAFF,  MS-F665 

PACIFIC  NORTHWEST  NATIONAL 
LABORATORY 
P.  O.  BOX  999 
BATTELLE  BOULEVARD 
RICHLAND,  WA  99352 

ATTN:  D.  N.  HAGEDORN,  MS  K5-12 

SANDIA  NATIONAL  LABORATORIES 
ATTN:  MAIL  SERVICES 
P.  O.  BOX  5800 

ALBUQUERQUE,  NM  87185-0459 

ATTN:  TECHNICAL  STAFF,  DEPT  5704, 
MS  0655 

ATTN:  TECHNICAL  STAFF,  DEPT  5704, 
MS  0979 

ATTN:  TECHNICAL  STAFF,  DEPT  5736, 
MS  0655 

ATTN:  TECHNICAL  STAFF.  DEPT  931 1 . 
MS  1159 
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SEISMOLOGICAL  LABORATORY  252-21 
CALIFORNIA  INSTITUTE  OF  TECHNOLOGY 
PASADENA,  CA  81125 
ATTN:  T.  AHERNS 

UNIVERSITY  OF  CALIFORNIA 
EARTH  SCIENCE  DIVISION 
470  MCCONE  HALL,  LBNL,  90-2106 
BERKELEY.  CA  94720 

ATTN:  L.  JOHNSON,  MS-90-1 116 

DEPARTMENT  OF  THE  AIR  FORCE 


HEADQUARTERS 

AFTAC/TTR 

1 030  SOUTH  HIGHWAY  A1A 
PATRICK  AFB,  FL  32925-3002 
ATTN:  V.  HSU 


AIR  FORCE  RESEARCH  LABATORY 
29  RANDOLPH  ROAD 
HANSCOM  AFB,  MA  01731-5000 
ATTN:  DR  D.  HARKRIDER 
ATTN:  DR  D.  REITER 
ATTN:  GPE,  J.  LEWKOWICZ 
ATTN:  J.  RONG-SONG 

AIR  FORCE  RESEARCH  LABORATORY 
5  WRIGHT  STREET 
HANSCOM  AFB,  MA  01731-3004 
ATTN:  RESEARCH  LIBRARY,  TL 

AIR  FORCE  TECH  APPLICATIONS  CENTER 
1300  17TH  STREET,  SUITE  1450 
ARLINGTON,  VA  22209 
ATTN:  R.  BLANDFORD 

AIR  FORCE  TECHNICAL  APPLICATIONS 
CTR/TT 

1 030  S  HIGHWAY  AIA 
PATRICK  AFB,  FL  32925-3002 
ATTN:  CA/STINFO 
ATTN:  DR  B.  KEMERAIT 
ATTN:  DR  D.  RUSSELL 
ATTN:  G.ROTHE.TTR 
ATTN:  J.C.  LUCAS 
ATTN:  M.SIBOL 

DIRECTORY  OF  OTHER  (LIBRARIES  AND 
UNIVERSITIES) 

ARIZONA,  UNIVERSITY  OF 
DEPT.  OF  G  EOSCI ENCES/SASO 
TUCSON,  A Z  85721 

ATTN:  PROF  T.  C.  WALLACE 


BOISE  STATE  UNIVERSITY 
GEOSCIENCES  DEPARTMENT 
1910  UNIVERSITY  DRIVE 
BOISE,  ID  83725 

ATTN:  J.  E.  ZOLLWEG 


BOSTON  COLLEGE 
INSTITUTE  FOR  SPACE  RESEARCH 
140  COMMONWEALTH  AVENUE 
CHESTNUT  HILL,  MA  02167 
ATTN:  PROF  L.  SYKES 


CALIFORNIA  INSTITUTE  OF  TECHNOLOGY 
DIVISION  OF  GEOLOGY  &  PLANETARY 
SCIENCES 

PASADENA,  CA  01125 

ATTN:  PROF  D.  V.  HELMBERGER 


CAUFORNIA-DAVIS,  UNIVERSITY  OF 
DAVIS,  CA  95616 

ATTN:  R.  H.  SHUMWAY,  DIV  STATISTICS 
ATTN:  S.E.  KIM 


CALIFORNIA-SANTA  CRUZ,  UNIVERSITY  OF 
INSTITUTE  OF  TECTONICS 
SANTA  CRUZ,  CA  95064 
ATTN:  DRR.  S.WU 
ATTN:  PRO  T.  LAY 

COLORADO-BOULDER,  UNIVERSITY  OF 
BOULDER,  CO  80309 

ATTN:  M.  RITZWOLLER, 

CAMPUS  BOX  390 
ATTN:  PROF  C.  ARCHAMBEAU 

COLUMBIA  UNIVERSITY 
LAMONT-DOHERTY  EARTH  OBSERVATORY 
PALISADES,  NY  10964 
ATTN:  DR  L.  R.  SYKES 
ATTN:  DRJ.XIE 
ATTN:  PROF  P.  G.  RICHARDS 

CONNECTICUT,  UNIVERSITY  OF 
DEPT.  OF  GEOLOGY  &  GEOPHICS 
STOORS,  CT  06269-2045 

ATTN:  PROF  V.  F.  CORMIER,  U-45, 
ROOM  207 


CORNELL  UNIVERSITY 

DEPARTMENT  OF  GEOLOGICAL  SCIENCES  . 

3126  SNEE  HALL 

ITHACA,  NY  14853 

ATTN:  PROF  M.  BARAZANGI 
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HARVARD  UNIVERSITY 
HOFFMAN  LABORATORY 
20  OXFORD  STREET 
CAMBRIDGE,  MA  02138 

ATTN:  PROF  A.  D2IEWONSKI 
ATTN:  PROF  G.  EKSTROM 


IRIS 

1200  NEW  YORK  AVENUE,  NW,  SUITE  800 
WASHINGTON,  DC  20005 
ATTN:  DR  D.  SIMPSON 
ATTN:  DR  G.  E.  VAN  DER  VINK 

MASSACHUSETTS  INSTITUTE 
OF  TECHNOLOGY 
EARTH  RESOURCES  LABORATORY 
42  CARLETON  STREET 
CAMBRIDGE,  MA  02142 

ATTN:  PROF  M.  N.  TOKSOZ 

MICHIGAN  STATE  UNIVERSITY  LIBRARY 
450  ADMINISTRATION  BUILDING 
EAST  LANSING,  Ml  48824 
ATTN:  K.  FUJITA 

NEW  MEXICO  STATE  UNIVERSITY 
DEPARTMENT  OF  PHYSICS 
LAS  CRUCES,  NM  88003 
ATTN:  PROFJ.NI 
ATTN:  PROF  T.  HEARN 

PENNSYLVANIA  STATE  UNIVERSITY 
GEOSCIENCES  DEPARTMENT 
403  DEIKE  BUILDING 
UNIVERSITY  PARK,  PA  16802 
ATTN:  PROF  C.  A.  LANGSTON 
ATTN:  PROFS. ALEXANDER 

SAN  DIEGO  STATE  UNIVERSITY 
DEPT  OF  GEOLOGICAL  SCIENCES 
SAN  DIEGO,  CA  92182 
ATTN:  PROF  S.  M.  DAY 

SOUTHERN  METHODIST  UNIVERSITY 
FONDREN  LIBRARY 
DALLAS,  TX  75275 
ATTN:  B.  STUMP 

ATTN:  G.  MCCARTOR,  DEPT  OF  PHYSICS 
ATTN:  H.  L.  GRAY,  DEPT  OF  STATISTICS 

UNIVERSITY  OF  COLORADO 
CAMPUS  BOX  583 
BOULDER,  CO  80309 

ATTN:  DR  A.  L.  LEVSHIN 


UNIVERSITY  OF  SOUTHERN  CALIFORNIA 
520  SEAVER  SCIENCE  CENTER 
UNIVERSITY  PARK 
LOS  ANGELES,  CA  90089-0483 
ATTN:  PROF  C.  G.  SAMMIS 


FOREIGN 


AUSTRALIAN  GEOLOGICAL 
SURVEY  ORGANIZATION 
CORNER  OF  JERRAGOMRRRA  AND 
NINDMARSH  DRIVE 
CANBERRA,  ACT  2609 
AUSTRALIA 

ATTN:  D.JEPSON 

GEOPHYSICAL  INSTITUTE  OF  ISRAEL 
HAMASHBIR  STREET,  1 
HOLON,  58122  ISRAEL 

ATTN:  DR  Y.  GITTERMAN 

I.R.I.G.M.  -  B.P.  68 
38402  ST.  MARTIN  D'HERES 
CEDEX,  FRANCE 

ATTN:  DR  M.  BOUCHON 

MINISTRY  OF  DEFENSE/ 
PROCUREMENT  EXECUTIVE 
ALACKNESS,  BRIMPTON 
READING  FG7-4RS  ENGLAND 
ATTN:  DR  P.  MARSHALL 

NTNF/NORSAR 
P.O.  BOX  51 

N-2007  KJELLER.  NORWAY 
ATTN:  DR  F.  RINGDAL 
ATTN:  T.  KVAERNA 


OBSERVATOIRE  DE  GRENOBLE 
I.R.I.G.M.  -  B.P.  53 
38041  GRENOBLE,  FRANCE 
ATTN:  DR  M.  CAMPILLO 

RESEARCH  SCHOOL  OF  EARTH  SCIENCES 
INSTITUTE  OF  ADVANCES  STUDIES 
G.  P.  O.  BOX  4 

CANABERRA  2601 ,  AUSTRALIA 
ATTN:  PROF  B.  L.  N.  KENNETT 

RUHR  UNIVERSTY/BOCHUM 
INSTITUTE  FOR  GEOPHYSIK 
P.  O.  BOX  102148 
463  BOCHUM  1,  GERMANY 
ATTN:  PROF  H.  P.  HARJES 
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SEISMOLOGICAL  DIVISION 
IRPG 

P.O.BOX  2280 
HOLON  58122 
ISRAEL 

ATTN:  A.  SHAPIRA 

SOCIETE  RADIOMANA 
27  RU  CLAUDE  BERNARD 
75005  PARIS.  FRANCE 
ATTN:  DR  B.  MASSINON 
ATTN:  DR  P.  MECHLER 

UNIVERSITY  OF  BERGEN 
INSTITUTE  FOR  SOUD  EARTH  PHYSICS 
ALLEGATION  40 
N-5007  BERGEN,  NORWAY 
ATTN:  R.E.HUSEBYE 

UNIVERSITY  OF  CAMBRIDGE 
DEPARTMENT  OF  EARTH  SCIENCES 
MADINGLEY  RISE,  MADINGLEY  ROAD 
CAMBRIDGE  CB3  OEZ,  ENGLAND 
ATTN:  PROF  K.  PRIESTLEY 

OTHER  GOVERNMENT 

CENTRAL  INTELLIGENCE  AGENCY 
WASHINGTON.  DC  20505 

ATTN:  CHIEF.  OTI/MNG/NWTT  -  ST47  NHB 

NATIONAL  ARCHIVES  AND 
RECORDS  ADMINISTRATION 
8601  ADELPHI  ROAD,  ROOM  3360 
COLLEGE  PARK,  MD  20740-6001 
ATTN:  USER  SERVICE  BRANCH 

NATIONAL  PHOTOGRAPHIC 
INTERPRETATION  CENTER 
SOUTHWEST  STATION 
P.  O.  BOX  70967 
WASHINGTON,  D  C  20024-0967 
ATTN:  K.  TIGHE-WHITE 
ATTN:  M.  BURNS 

U  S  DEPARTMENT  OF  STATE 
320  21  ST  STREET,  NW 
WASHINGTON,  DC  20451 

ATTN:  DR  E.  LACEY,  ACDA/IVI, 

ROOM  5741 
ATTN:  K.  WARD 
ATTN:  R.  DAY 
ATTN:  T.  RAY 


US  DEPARTMENT  OF  THE  INTERIOR 
US  GEOLOGICAL  SURVEY  NATIONAL  CENTEI 
MILITARY  GEOLOGY  PROJECT 
12201  SUNRISE  VALLEY  DRIVE 
RESTON,  VA  22029 
ATTN:  B.  LEITH 
ATTN:  DRJ.FILSON 
ATTN:  W.  LEITH,  MS  928 
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CUT  HERE  AND  RETURN 


DISTRIBUTION  LIST  UPDATE 


This  mailer  is  provided  to  enable  DTRA  to  maintain  current  distribution  lists  for  reports.  fWe  would 
appreciate  vou  providing  the  requested  information^ 

□  Add  the  individual  listed  to  your  distribution  list 

□  Delete  the  cited  organization/individual. 

□  Change  of  address. 

NAME: _ 

'ORGANIZATION: _ _ _ 

OLD  ADDRESS  NEW  ADDRESS 


Note: 

Please  return  the  mailing  label  from  the 
document  so  that  any  additions,  changes, 
corrections  or  deletions  can  be  made  easily. 
For  distribution  cancellation  or  more 
information  call  DTRA/ADM  (703)  767-4724. 


TELEPHONE  NUMBER:  (  ) _ _ 

DTRA  PUBLICATION  NUMBER/TITLE  CHANGES/DELETIONS/ADDITIONS,  etc.) 

(Attach  Shaet  if  more  Spaca  is  Raqulnd) 


DTRA  or  other  GOVERNMENT  CONTRACT  NUMBER: _ 

CERTIFICATION  of  NEED-TO-KNOW  BY  GOVERNMENT  SPONSOR  (if  other  than  DTRA): 
SPONSORING  ORGANIZATION: 

CONTRACTING  OFFICER  or  REPRESENTATIVE: 

SIGNATURE: 


DEFENSE  THREAT  REDUCTION  AGENCY 
ATTN:  ADM 

8725  John  J.  Kingman  Road,  MS-6201 
Ft  Belvoir,  VA  22060-6201 


DEFENSE  THREAT  REDUCTION  AGENCY 
ATTN:  ADM 

8725  John  J.  Kingman  Road,  MS-6201 
Ft  Belvoir.  VA  22060-6201 


